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Characterization of qubit couplings in many-body quantum systems is essential for benchmarking 
quantum computation and simulation. We propose a tomographic measurement scheme to determine 
all the coupling terms in a general many-body Hamiltonian with arbitrary long-range interactions, 
provided the energy density of the Hamiltonian remains finite. Different from quantum process 
tomography, our scheme is fully scalable with the number of qubits as the required rounds of 
measurements increase only linearly with the number of coupling terms in the Hamiltonian. The 
scheme makes use of synchronized dynamical decoupling pulses to simplify the many-body dynamics 
so that the unknown parameters in the Hamiltonian can be retrieved one by one. We simulate the 
performance of the scheme under the influence of various pulse errors and show that it is robust to 
typical noise and experimental imperfections. 


Introduction. —Physicists have been striving to under¬ 
stand and harness the power of quantumness since the 
establishment of the quantum theory. With the flour¬ 
ishing of quantum information science in recent decades 
[1, 2], numerous breakthroughs—both in theory and 
in experiment—helped to frame a clearer goal: it is 
the entanglement and the exponentially growing Hilbert 
space that distinguishes quantum many-body systems 
from classical systems [3-5]. To fully leverage the quan¬ 
tum supremacy, a vital step is to verify and benchmark 
the quantum device. The standard techniques of quan¬ 
tum state and process tomography [6-11], however, are 
plagued by the same exponential growth of dimensions 
[12]. A related problem is to directly identify Hamilto¬ 
nians, the generators of quantum dynamics. They can 
often be specified by fewer number of parameters that 
scales polynomially with the system size. 

Hamiltonian tomography for generic many-body sys¬ 
tems is nevertheless a daunting task. The way to extract 
information of unknown parameters in a Hamiltonian is 
by measuring certain features of its generated dynam¬ 
ics. To make this possible, one has to solve the dynamics 
generated by the Hamiltonian to make a definite connec¬ 
tion between its dynamicai features and the Hamiitonian 
parameters. However, for generai many-body Hamii- 
tonians, their dynamics are extremeiy compiicated and 
intractabie by numericai simuiation as the simuiation 
time increases exponentiaiiy with the size of the system. 
Progress in this direction has mostiy be on smaii sys¬ 
tems [13-17] or speciai many-body systems which are ei¬ 
ther exactiy soivabie due to many conserved operators, of 
limited Hilbert space dimensions amenable to numerical 
simulation, or short-range interacting systems [18-24]. 

In this paper, we propose a scheme to achieve Hamil¬ 
tonian tomography for general many-body Hamiltonians 
with arbitrary long-range couplings between the qubits. 
The key idea is to simplify the dynamics generated by a 
general many-body Hamiltonian through application of 


a sequence of dynamical decoupling pulses on individual 
qubits. Dynamical decoupling (DD) is a powerful tech¬ 
nique that uses periodic fast pulses to suppress noise and 
average out unwanted couplings between the system and 
the environment [25-40]. We apply a sequence of syn¬ 
chronized DD pulses on a pair of qubits, which forms a 
small target system that has coupling with the rest of the 
qubits in the many-body Hamiltonian, the effective en¬ 
vironment. The DD pulses keep the desired couplings 
within this target system intact while average out its 
couplings with all the environment qubits. The dynam¬ 
ics under the DD pulses become exactly solvable, from 
which we can perform a tomographic measurement to de¬ 
termine the coupling parameters within this small target 
system [13-16]. We then scan the DD pulses to different 
pairs of qubits to measure all the other coupling terms 
in the Hamiltonian. We assume the ability to address 
individual qubits, which is realistic for many experimen¬ 
tal platforms, such as trapped ions [41-43], cold atoms 
[44-46], and solid-state qubit systems [47-50]. Several 
features make the scheme amenable to experimental im¬ 
plementation. First of all, applying the DD pulse se¬ 
quence is a standard procedure in many experiments. 
Post-processing of data is straightforward as it only re¬ 
quires one or two parameter curve fitting. In addition, 
we demonstrate with explicit numerical simulation that 
the scheme is robust to various sources of errors in prac¬ 
tical implementation, such as the remnant DD coupling 
error, measurement uncertainties, and different types of 
pulse errors. 

Scheme for Hamiltonian tomography. —The system we 
have in mind is the most general Hamiltonian with two- 
body qubit interactions 

H= E + (1) 

a,0 ,m<n m,a 

where characterizes the coupling strength between 
spins m and n for the a, (3 components, and 6^ represents 


2 





Time 

Y X Y Y X Y 



FIG. 1. (color online). Schematics for the tomography pro¬ 
cedure. (a) To map out the coupling coefficients J^f ^ a syn¬ 
chronized DD sequence is applied to spins i and j. Both spins 
will be decoupled from the rest of the system, (b) The XY- 
8 DD sequence on spins i and j to probe the parameters of 
the Hamiltonian in Eq. (3). The initial state is for instance 
prepared to the |00) state for the two spins, (c) To retrieve 
information about the local fields XY-8 pulse sequences 
are applied to the environment spins to decouple spin i from 
the rest. 


the local field on spin m; (a^) are the Pauli matrices 
along the a {/3) direction with a,/3 e {x^y^z). To adopt 
consistent notations throughout the text, we use m, n to 
denote a general spin label and i, j to refer to the specific 
target spins that we are probing with the DD pulses, 
calling the rest of the spins as environment spins. The 
terms spin and qubit are used interchangeably. Let the 
energy unit of the Hamiltonian be J, chosen to be the 
largest magnitude of all coefficients, so and h^/J 

are bounded between —1 and 1. In order to map out the 
coupling coefficient for the target spins, we propose 
to decouple these two spins from the environment spins 
by a synchronized DD pulse sequence. A synchronized 
XY-A sequence applied to both spins will average out 
their interactions with other spins while preserving the 
two-spin coherence (see Fig. l(a-b) for the schematic and 
the pulse sequence). Basically, only those interactions 
that commute with the DD sequence will survive. More 
rigorously, the evolution operator in one period is 

Ui = C/o'/VfcTjt/oafaJt/oafaJC/ocrfcTjf/o'/' 

where Uq = r is the time interval between two 

consecutive pulses, and 5, the bath, includes all terms of 
the Hamiltonian that only acts on environment spins. 
See Supplemental Material for the detailed derivation 
[51]. To bound the error term to O(J^r^), we assume 
~ interaction strength decays 

rapidly with spin separation distance so that the energy 
density of the Hamiltonian is bounded by a constant. 


This condition is satisfied for any finite systems as in the 
experiment with arbitrary interactions. In the thermo¬ 
dynamic limit, it is also a reasonable assumption for any 
physical systems whose energy is extensive. It may also 
be related to the generalized Lieb-Robinson bound for 
systems with long-range interactions [52-55]. The XY-8 
pulse sequence, which is the concatenation of XY-A se¬ 
quence with its time-reversal, eliminates the error term 
to the third order Fig. 1(b) shows the XY-8 

DD pulse sequence. Hence, in the Hilbert subspace of 
the target spins, the effective Hamiltonian is 

-ff 2 -spin = Cicrf ctJ + C 2 CrfcrJ + c^a-aj, (3) 

where we use Ci = jy, C 2 = JfP C 3 = jy to simplify the 
notation. The effective two-spin unitary evolution after 
Nc cycles of XY-8 sequence is 


t^ 2 -spin — 

cos((ci-C 2)T) 
gic3 T 

0 


0 

sin((ci-C2)T) 


0 

cos((ci+C2)T) 

e-icsT 

sin((ci+C2)T) 

ie — ic^T 

0 


0 

sin((ci+C2)T) 

^g-icgT 

cos((ci+C2)T) 
Q — ic^ T 

0 


sin((ci-C2)T) 

ieic^T 

0 

0 

cos((ci-C2)T) 

eic^T 


\ 
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where T = 8NcT is the total time. From the above ex¬ 
pression, one may notice that the Hamiltonian param¬ 
eters can be retrieved by preparing a particular initial 
state and measuring its time-evolved output probability 
in a given basis. In particular, we have 

P|+I)^| 00 > = |(00|[/2-spin|+I)|' = i [1 + sm(2(ci - C 2 )T)] 
-P|+I)-!-|10) = KlO|P2-spin|+I)|^ = i sin(2(ci + C 2 )r)] 
-P|0I)-!-|++) = K + +|P2-spin|0I)|^ = J [1 + sin(2(c2 — C 3 )r)] 


where |+) = ^(| 0 ) + | 1 )) and |I) = ^(| 0 ) +i|l)) are 
the rotated basis. The coupling strengths ci,C 2 and C 3 
can be extracted from the oscillation frequencies of these 
three sets of measurements at various time points. These 
particular sets are not the only suite to extract those 
parameters. They are chosen for the convenience in fit¬ 
ting and in state preparation. Only product states of the 
two target spins, disentangled from the rest, are required. 
We also remark that the error incurred is 0{NcJ^r^), so 
one needs Jr 1 for a robust decoupling scheme. In a 
similar fashion, one can retrieve all other coupling coef¬ 
ficients. Let us denote the synchronized XY-8 DD pulse 
sequence as XiXj-YiYj-8 to show explicitly the particu¬ 
lar pulses on specific spins. Replacing the sequence with 
XiYj-YiZj-8 (YiXj-ZiYj-8) pulses, we will be able to ex¬ 
tract the coefficients , jy^ and ^ Jif and 

respectively. 

By scanning the DD pulses to different target pairs, 
the above procedure recovers all the coupling coefficients 
retrieval of local field coefficients follows a 
similar approach. We now need to decouple the par¬ 
ticular spin i from the rest without contaminating its 
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FIG. 2 . (color online). Numerical simulation and curving fitting results, (a)-(c) are used to retrieve , Jjg and J 7 I between 
spins 7 and 9. (d) and (e) are used to extract and bl for spin 6 . Each measurement data prn are drawn from the binomial 

distribution with the true probability p as the mean and p(l — p)/Nm as the variance. The measurement uncertainty of each 
point is thus y^prn(l — Pm)/Nm- The blue solid lines are the best-fit lines with the simulated experimental data prn, and red 
dashed lines are the theoretical ones generated by the true Hamiltonian parameters. Pulse errors are not included in these 
plots, so any discrepancies stem from the remnant DD coupling error and measurement uncertainties. Other parameters used 
are N = 12, rJ = 0.01, Nm = 100, Nt = 50. 


own spin term Shining a XF -8 DD sequence 

on spin i removes all information about 69 too. In¬ 
stead, one could address all the environment spins with 
XY-^ pulses, and decouple them from spin i (alterna¬ 
tive schemes are discussed in the Supplemental Material 
[51]). This scheme will be very robust to pulse errors, 
since no laser pulses are directly applied to the target 
spin [Fig. 1(c)]. The effective single-spin Hamiltonian 
is thus ii7i-spin = with a unitary 

evolution 7/i-spin = g-*^i-spinT executing a spin rota¬ 
tion on the Bloch sphere. Again, by preparing a par¬ 
ticular state and measuring its time evolution, we get 


| 0 )^| 0 > 


1 + 


iPt/bf 


sin^( 6 T) and P|+)^|+) = 1 + 


{bf /b)‘^ — 1 sin^( 6 T), where 6 = \/ 


is the magnitude of the Bloch vector. These two sets of 
measurements will determine bf^b^ and 6 f up to a sign. 
The correct signs from the remaining discrete set can be 
picked out by measuring P|+)^|o) and P|i)^|o) at a single 
time point [51]. 


The complete scheme applies to any generic Hamilto¬ 
nian with interacting qubits. In the most general case, 
one needs to determine 9N{N — l)/2 + 3N coefficients. 
However, in many physical systems, the particular form 
of the interaction is known and/or the interaction of¬ 
ten decays fast enough that one can significantly reduce 
the number of measurements required. In particular, if 
can be truncated at some spin separation distance 
in the case of short-range interactions, the number of 
measurements will be linear with the system size N. In 
the following, we numerically simulate the experimental 
procedure for the most general Hamiltonian, taking into 
account various sources of errors, including the remnant 
DD coupling error, measurement uncertainties, and dif¬ 
ferent forms of pulse errors. 


Numerical simulation .—We consider the general 
Hamiltonian given in Eq. ( 1 ) with coefficients / J and 
69/J randomly drawn from —1 to 1. In our finite-system 


simulation, we ignore the decay of with distance, 
so the system may include unphysically long-range in¬ 
teractions and could simulate Hamiltonians in any di¬ 
mensions. To retrieve J9T, for example, we start with 
a product state of all spins, and perform time evolution 
using the entire Hamiltonian from Eq. (1), interspersed 
with the XY-3 DD pulses on target spins i and j. We 
would like to emphasize that specific state initialization 
for the environment spins is not required as long as they 
are disentangled from the target pair of qubits at the 
beginning. After Nc cycles of the DD sequence, the en¬ 
vironment spins are traced out and measurements are 
made on spins i and j. In the simulation, we do not as¬ 
sume the pure unitary evolution 6 / 2 -spin as the remnant 
coupling to the environment spins may entangle the two 
spins with the rest. However, any undesired couplings 
are suppressed to the order of O(J^r^) and we do ob¬ 
serve that the two-spin density matrix remains mostly 
pure (^ 99.9%) for our chosen parameters. 

As the tomography procedure involves measuring the 
output probability of a certain state, each time point 
will be measured Nm times, which gives an estimate of 
the probability Pm in this state. The measurement un¬ 
certainty (standard deviation) will be y^Pm{^ —Pm)/Nm 
following the binomial distribution. As discussed above, 
to map out Cl, C 2 and C 3 , one needs to measure 
^ 1 + 1 )^ 100 ), ^ 1 + 1 )^ 110 ) and P|oi)^|++) for the target spins 
at various time points and extract the corresponding os¬ 
cillation frequencies. Suppose Nf different time points 
are measured for each set. The oscillation frequencies 
can be found either by Fourier transform or by curve 
fitting. In general, if data show numerous oscillation 
periods, Fourier transform will be more robust and re¬ 
liable [14-16]. In our case, however, the long time ob¬ 
servations will be undermined by the remnant coupling 
to the environment spins and possible pulse error accu¬ 
mulation. Simple curving fitting with fewer oscillation 
periods, therefore, appears to be a better solution. In 
Fig. 2(a-c), we fit the data with the method of least 
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squares with rJ = 0.01, = 100, A^t = 50 for spins 

i = 7 and jf=9inaA^ = 12 spin system. The blue solid 
lines are the best-fit lines, and the red dashed lines are 
the theoretical lines using the true coupling coefficients. 
The longest time period requires 800 pulses, which is well 
within the current experimental technology without sig¬ 
nificant pulse error accumulation. Table I compares the 
true values and the estimated ones of . Uncertainties 
in the estimation stem from the curve fitting due to mea¬ 
surement uncertainties. Corresponding results for of 
spin 6 are shown in Fig. 2(d-e) and Table I. All estimated 
parameters are accurate within a few percent. 

To simulate real experiments, one also needs to include 
possible pulse errors. One possible source of errors is 
the finite duration of each control pulse, which limits 
the minimum cycle time. This is typically not the dom¬ 
inant source of errors and can often be well-controlled 
[31, 35, 56-58]. In most experiments, the major cause 
of errors is the deviation between the control pulses and 
the ideal X or Y pulses. These can either arise from the 
amplitude error where the rotation angle differs from the 
ideal 7r-pulse or the rotation error where the rotation axis 
deviates from the x or y axis. In typical experiments, in¬ 
dividual pulse errors may be controlled within a percent 
level. In our simulation, we consider three different forms 
of pulse errors: Systematic Amplitude pulse Error (SAE), 
Random Amplitude pulse Error (RAE) and Random Ro¬ 
tation axis Error (RRE). See the caption of Table I for 
the specific forms of the errors. Moderate systematic er¬ 
rors can be self-compensated by the XY-S DD sequence. 
Numerically, we found that 5% of SAE has negligible ef¬ 
fect on the parameter estimation. In addition, we also 
simulated the cases where each pulse experiences a 1% 
RAE or RRE. Results are summarized in Table I. The 
average deviation from the true parameters are within 
5%. Here, we would like to point out a few features of 
our scheme that make it inherently robust to errors. Eirst 
of all, the estimation of the coupling strength only 
entails frequency estimation, which could endure large 
deviations of a few measurement points. In addition, the 
single-parameter curve fitting scheme not only makes the 
estimation robust but is also more convenient for experi¬ 
ment. Moreover, the retrieval of local fields is remark¬ 
ably tolerant to pulse errors. Since no pulse is directly 
applied to the target spin, any pulse errors on the envi¬ 
ronment spins will only be propagated via the remnant 
DD coupling error, which is suppressed to the order of 
We have numerically tested that a 10% pulse 
error of any kind would have negligible effects on the es¬ 
timation of b^. Alternative schemes to extract the local 
fields are detailed and discussed in the Supplemental Ma¬ 
terial [51]. They are less tolerant to pulse errors, but may 
be easier to implement in some experimental setups. 

Discussion and outlook .—We have thus numerically 
demonstrated that the proposed scheme is robust to var¬ 
ious sources of errors present in real experiments. The 


TABLE I. NPE: No Pulse Error; SAE: Systematic Amplitude 
pulse Error; RAE: Random Amplitude pulse Error; RRE: 
Random Rotation axis Error; AD: Average Deviation from 
true values. The last digit in bracket for each number quan¬ 
tifies the estimation error bar due to measurement uncertain¬ 
ties, which is generated by the bootstrapping method. The 
percentage values in the brackets denote the amount of errors 
introduced in each pulse. The errors are in the form of: SAE, 

where e = 5%, 6 is randomly chosen from (—1%, 1%), (<a, 7 ) 

is a vector with a random direction but fixed magnitude at 
1%, and n — x^y for the X and Y pulses respectively. 



True 


Estimated Parameters 



- 

NPE 

SAE(5%) 

RAE( 1 %) 

RRE(1%) 

JXX 

^79 

-0.378 

-0.369(3) 

-0.377(3) 

-0.379(4) 

-0.412(2) 

jyy 

^79 

0.863 

0.856(3) 

0.846(3) 

0.867(4) 

0.836(2) 

JZZ 

^79 

0.679 

0.669(5) 

0.649(5) 

0.718(6) 

0.611(4) 

h% 

0.334 

0.32(1) 

0.32(1) 

0.32(1) 

0.32(1) 

K 

0.569 

0.567(8) 

0.567(8) 

0.567(8) 

0.568(8) 

hi 

-0.431 

-0.441(8) 

-0.443(8) 

-0.441(8) 

-0.441(8) 

AD 

- 

2 % 

3% 

3% 

5% 


measurement uncertainties can be lowered by increasing 
Nm and the pulse errors can be reduced by limiting the 
maximum number of pulses needed. The optimal strat¬ 
egy involves a delicate balance between experimental so¬ 
phistication and error control. Eor example, by fixing 
tJ and the total number of measurements for each set, 
Nm X A^t, one could devise an optimal estimation proce¬ 
dure. In addition, it is also possible to eliminate the rem¬ 
nant DD coupling error to a higher order with more elab¬ 
orate pulse sequences such as the concatenated DD se¬ 
quence [30, 31] and reduce pulse errors by designing com¬ 
posite pulses or self-correcting sequences [34, 35, 59, 60]. 
The scheme can also be extended straightforwardly to 
qudit systems of higher spins or to bosonic or fermionic 
systems. 

In conclusion, we have proposed a general scheme to 
achieve full Hamiltonian tomography for generic interact¬ 
ing qubit systems with arbitrary long-range couplings. 
The required number of measurements scales linearly 
with the number of terms in the Hamiltonian, and the 
scheme is robust to typical experimental errors or imper¬ 
fections. 
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In this Supplemental Material, we provide more details on the dynamical decoupling scheme and error 
estimation, the retrieval of local field parameters and alternative schemes, and include further results 
taking into account of pulse errors. 


DYNAMICAL DECOUPLING 


The most general Hamiltonian with two-body qubit interactions can be written as 


a,/3,m<n t 


|Ck: a 

m • 


(4) 


The energy unit of the Hamiltonian is taken to be J such that and h^/ J are bounded between —1 and 1. The 

symmetric XY-A dynamical decoupling (DD) sequence on both spins i and j produces 


Ui = (crfcrjt/ocrfaj) «(7|C/ocrf {afa^Uoafa^j) (5) 

where Uq = and r is the time interval between consecutive pulses. We can decompose H into two parts. 


rl/2 


^ = ^0 + ^^!, 

Hi = 

+ uf + afBf + alBt + aJS* + a^B^ + a]B^ + B, 


(6) 

(7) 

(8) 


where Bf includes the local field on the rth spin and interacting terms between the rth spin and all other spins other 
than the jth spin, i.e., Bf = hf + n^i j •^iXn- bath term B includes all the environment operations, i.e., all 
operators that does not act on spins i and j. We define other Hamiltonian part as 


afafHafa] = Ho + H2, ufajH-ufuf = Fq + 7^3, = Ho + H4, 


* 3 


(9) 


where 


' ^3 ^3 




' ^3 V 


'ij ^j 


+ uf Bf 


af Bf - af Bf + aJBJ - ujBj - uf Bf + B, 
Hi = -J-fafa] + Jifcrfa] - - -^f/ufcrf + 

- erf Bf + afBf - ufBf - uf Bf + erf Bf - erf Bf + B, 
H. = +JfFafa'i - JYafY + - JffafY - JYo 


ZX 

ij 


* J 


. J^yrr^rry 




- afBf - afBf + a^B^ - a]B- - a^^B^ + uf Bf + B. 


( 10 ) 

( 11 ) 

( 12 ) 


Basically, each term will either commute or anticommute with the operator afa^. Those commuting with it will be 
left invariant, and those anticommuting will have a flipped sign. Hq and B commute with each operator so 

they are left unchanged. Now we can see explicitly that Hi Y H 2 Y Y = 4H, which is why the DD sequence 
effectively decouples the two spins i and j with the rest of the spins. To estimate the error, we combine the unitary 
evolution for a period and repeatedly make use of the formula 


(13) 


Ignoring and higher-order terms, we find 


= ^-ir/ 2 (Ho+Hi) ^-iT(Ho+H2) ^-iT/ 2 (Ho + Hi) 

^ ^-iAT(Ho+B)+C 


(14) 



where the remnant coupling noise term is 


C = -\t^ [Ho + Fi, ffo + H 2 ] - [|ffo + + ^^2, Ho + H^] 

- P" [iHo + \Hi +H2+ Hi, Ho + Ho] - y [iHo + \Hi +H2 + H0 + H^, Ho + H^] + 0{t^) 

= - t2 {[^0,^3 - H 2 ] + \[H2 - Ho,Hi] + \[H2,Ho\}+0{t^). (15) 

In the error term C, the biggest contribution comes from terms like Our aim is to show that the error 

does not scale with the system size i.e., C = Let us consider one such term and write it out explicitly 

(suppressing the a^jS summation): 




E ja^ a B ^ 

i / j ^ip p 

m<n PT^iJ 


m<n 


(16) 


Since and Jfj are rapidly decaying functions of the separation distance, for a fixed site X]m<n ^rnnJim ~ 

Note that this differs from the scaling of the Hamiltonian, H ^ Jmn — 0{N J). All the other terms in C are 

either smaller or contribute to the same order as the above term. Therefore, we have C = 0{J‘^r‘^). To be able to 
neglect the error terms, one needs to fulfill the condition Jr <^1. 

The above discussion is pertinent to the XY-A pulse sequence. We can cancel the second order contribution by 
using the XT-8 pulse sequence as U 2 = UiUf^^ where is just the time-reversed sequence of Ui. It can be readily 
seen that the error terms C and will cancel each other to the second order O(r^), since contains the same 
terms as in C only with the role of H 2 and Hs interchanged. Therefore, the remnant coupling error of the XT-8 pulse 
sequence is as discussed in the main text. 


LOCAL FIELD RETRIEVAL 


In the main text, we proposed a scheme to retrieve the local fields bf by shining the XT-8 pulse sequences on all 
the environment spins. Here, we provide more details and outline alternative schemes that may in some experimental 
setups be easier to implement. By decoupling the environment spins with spin i as illustrated in Fig. 1(c) of the main 
text, we have the effective single-spin Hamiltonian iLi-spin = bfaf + + b^af. The time evolution operator is 


Uu 


spin 




/ b^ ib^ ^by \ 

cos(6 T)-z^sin(6 T) - ^ ^ ^ sin(6r) ' 

bj — ib^ . , 6 : 


(17) 


——- sin(6T) cos(6T) + sin(6T) j 

where b = \/{b^P + {b\Y + (6|)^ is the magnitude of the Bloch vector. By measuring 


| 0 )^| 0 ) 


l+>^l+> 


= 1 


= 1 + 


{bl/bf-l sin2(6T) 
{b^bf - ll sin2(6T) 


(18) 

(19) 


at various time points, we could determine \bf\, \b^\, |6||. To pin down the correct signs, one can supplement the above 
two sets of measurements with another two measurement points; 


^+)^|0) — |(0|17l-spin|+)|^ “ 9 ( ^ 




P|i)^|o) = |(0|f/i..pi„|l)|" = i(l + 


&2 

2bU: 




i2bT 


ip-siiP bT+^sm2bT 
b^ b 


( 20 ) 

( 21 ) 


Only one time point is needed to determine the signs. For example, one could take measurements at bT = 7r/4 and 
use P|+)^|o) and P|q^|o) to pick out the correct signs. 

The above procedure requires applying the DD sequences to all spins other than the target spin. In some experi¬ 
mental setting, it may be easier to apply a global DD sequence to all spins and add another individually addressed 
beam on spin i to cancel the DD sequence on that single spin. See Fig. 3 for illustration. For instance, one could apply 
synchronized XahIah-S global pulses and in addition XiYi-8 focused pulses on spin i. In this way, spin i effectively 
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m 4 0 0 ^ 

FIG. 3. Alternative scheme to map out the local fields b^. A global pulse imposes the XY-8 pulse sequence on all spins and a 
focused pulse is in addition applied to spin i to cancel the DD sequence on that single spin. 


experiences no pulses at all time. The effective Hamiltonian again reduces to the same i^i-spin as above. However, 
this scheme is not very robust to pulse errors. Any deviation from the ideal pulse will be doubled on spin i and 
accumulate. The pulse error will affect the single-spin coherence and obscure bf too. We have tested it numerically 
that the pulse errors have to be controlled within 0.5% for the scheme to be feasible. So it can be used in some setups 
where pulse errors are not an issue or the total number of pulses can be reduced. One may also use this scheme and 
modify the sequence by designing composite pulses or self-correcting sequences to reduce pulse errors. 


PULSE ERRORS 


In the main text, we discussed different types of pulse errors. In our numerical simulation, we considered Systematic 
Amplitude pulse Error (SAE), Random Amplitude pulse Error (RAE) and Random Rotation axis Error (RRE). The 
fitting curves in Eig. 2 of the main text do not take into account of pulse errors. Here, we include the figures (Eig. 4) 
for the case with a 1% RRE. We can see, for example in Eig. 4(a), that the frequency estimation is still very accurate 
while some measurement points may have a notable mismatch. We may also notice that the estimation of bf is 
exceptionally robust to pulse errors since no pulse is applied to spin i in the scheme. Other pulse errors have similar 
effects on the estimation of parameters. 



FIG. 4. (color online). Numerical simulation and curving fitting results with a Random Rotation axis Error (RRE) for each 
pulse. The RRE is of the form +/3ay-\-ja ) is a vector with a random direction but fixed magnitude 

at 1%. (a)-(c) are used to retrieve J^q ^ J 79 and J 7 I between spins 7 and 9. (d) and (e) are used to extract ^6 foi* 

spin 6 . The blue solid lines are the best-ht lines with the simulated experimental data, and red dashed lines are the theoretical 
ones generated by the true Hamiltonian parameters. Other parameters are the same as in Fig. 2 of the main text. 







